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ABSTRACT 

We study the formation of dark matter halos in the concordance ACDM model over a wide range of redshifts, 
from z = 20 to the present. Our primary focus is the halo mass function, a key probe of cosmology. By performing 
a large suite of nested-box A^-body simulations with careful convergence and error controls (60 simulations with 
box sizes from 4 to 256/i~'Mpc), we determine the mass function and its evolution with excellent statistical 
and systematic errors, reaching a few percent over most of the considered redshift and mass range. Across the 
studied redshifts, the halo mass is probed over 6 orders of magnitude (10^ - lO'^^/z'^M©). Historically, there has 
been considerable variation in the high redshift mass function as obtained by different groups. We have made a 
concerted effort to identify and correct possible systematic errors in computing the mass function at high-redshift 
and to explain the discrepancies between some of the previous results. We discuss convergence criteria for the 
required force resolution, simulation box size, halo mass range, initial and final redshifts, and time stepping. 
Because of conservative cuts on the mass range probed by individual boxes, our results are relatively insensitive 
to simulation volume, the remaining sensitivity being consistent with extended Press-Schechter theory. Previously 
obtained mass function fits near z = 0, when scaled by linear theory, are in good agreement with our results at all 
redshifts, although a mild redshift dependence consistent with that found by Reed et al. may exist at low redshifts. 
Overall, our results are consistent with a "universal" form for the mass function at high redshifts. 

Subject headings: methods: A^-body simulations — cosmology: halo mass function 



L INTRODUCTION 

A broad suite of astrophysical and cosmological observations 
provides compelling evidence for the existence of dark matter 
Although its ultimate nature is unknown, the large-scale dy- 
namics of dark matter is essentially that of a self-gravitating 
collisionless fluid. In an expanding universe, gravitational in- 
stability leads to the formation and growth of structure in the 
dark matter distribution. The existence of localized, highly 
overdense dark matter clumps, or halos, is a key prediction of 
cosmological nonlinear gravitational collapse. The distribution 
of dark matter halo masses is termed the halo mass function 
and constitutes one of the most important probes of cosmol- 
ogy. At low redshifts, z < 2, the mass function at the high- 
mass end (cluster scales) is very sensitive to variations in cos- 
mological parameters, such as the matter content of the Uni- 
verse f^m, the dark energy content alon g with its equation-of- 
state parameter, w dHolder et al. Il200l1) . and the normalization 
of the primordial fluctuation power spectrum, ag. At higher 
redshifts, the halo mass functio n is important in probin g quasar 
abundance and formation sites (iHaiman & Loeb 11200 Ih. as well 
as the reionization history of the Universe ( Furlanetto et al. I 
120061) . 

Many recently suggested reionization scenarios are based on 
the assumption that the mass function is given reliably by mod- 
ified Press-Schechter type fits (Press & Schechter 1974, here- 
after PS; Bond et al. 1991). However, the theoretical basis of 
this approach is at best heuristic and careful numerical studies 
are required in order to obtain accurate results. Two examples 
serve to illustrate this statement. Reed et al. (2003) report a 
discrepancy with the Sheth-Tormen fit (Sheth & Tormen 1999, 
hereafter ST) of ~50% at a redshift of z = 15 (we explain the 



different fitting formulae and their origin in Heitmann et 
al. (2006a) show that the Press-Schechter form can be severely 
incorrect at high redshifts: at z > 10, the predicted mass func- 
tion sinks below the numerical results by an order of magnitude 
at the upper end of the relevant mass scale. Consequently, in- 
correct, or at best imprecise, predictions for the reionization 
history can result from the failure of fitting formulae. 

Since halo formation is a complicated nonlinear gravitational 
process, the current theoretical understanding of the mass, spa- 
tial distribution, and inner profiles of halos remains at a rela- 
tively crude level. Numerical simulations are therefore crucial 
as drivers of theoretical progress, having been instrumental in 
obtaining impo rtant results such as the Navarro-Frenk- White 
(NEW) profile (iNavarro et al. l[T997l) for dark matter halos and 
an (approximate) universal form for the mass function (Jenkins 
et al. 2001, hereafter Jenkins). In order to better understand 
the evolution of the mass function at high redshifts, a number 
of numerical studies have been carried out. High-redshift simu- 
lations, however, suffer from their own set of systematic issues, 
and simulation results can be at considerable variance with each 
other, differing on occasion by as much as an order of magni- 
tude! 

Motivated by all of these reasons we have carried out a nu- 
merical investigation of the evolution of the mass function with 
the aim of attaining good control over both statistical and, more 
importantly, possible systematic errors in A^-body simulations. 
Our first results have been reported in condensed form in Heit- 
mann et al. (2006a). Here we provide a more detailed and 
complete exposition of our work, including several new results. 

We first pay attention to simulation criteria for obtaining ac- 
curate mass functions with the aim of reducing systematic ef- 
fects. Our two most significant points are that simulations must 
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be started early enough to obtain accurate results and that the 
box sizes must be large enough to suppress finite-volume ar- 
tifacts. As in most recent work following that of Jenkins, we 
define halo masses using a friends-of-friends (FOF) halo finder 
with linking length b = 0.2. This choice introduces systematic 
issues of its own (e.g., connection to spherical overdensity mass 
as a function of redshift), which we touch on as relevant below. 
As it is not quantitatively significant in the context of this paper, 
we leave a detailed discussion to later w ork (Z. Lukic et al., in 
preparation; see also lReed et al. Il2007h . 

The more detailed results in this paper enable us to study the 
mass function at statistical and systematic accuracies reaching 
a few percent over most of our redshift range, a substantial im- 
provement over most previous work. At this level we find dis- 
crepancies with the "universal" fit of Jenkins at low redshifts 
(z < 5), but it must be kept in mind that the universality of the 
origin al fit was only meant to be at the ±20% level. Recently, 
iReed etal. I (120071) have found violation of universality at high 
redshifts (up to z = 30). To fit the mass function they have incor- 
porated an additional free parameter, the effective spectral in- 
dex Heft, with the aim of understanding and taking into account 
the extra redshift dependence missing from conventional mass- 
function-fitting formulae. Our simulation resu lts are consistent 
with the trends found by iReed et al. I (|2007|) at low redshifts 
(z < 5), but at higher redshifts we do not observe a statistically 
significant violation of the universal form of the mass function. 

Results from some previous simulations have reported good 
agreement with the Press-Schechter mass function at high red- 
shifts. Since the Press-Schechter fit has been found signifi- 
cantly discrepant with low-redshift results (z < 5), this would 
imply a strong disagreement with extending the well-validated 
low-redshift notion of (approximate) mass function universal- 
ity to high z. Our conclusion is that the simulations on which 
these findings were based violated one or more of the criteria to 
be discussed below. 

As simulations are perforce restricted to finite volumes, the 
obtained mass function clearly cannot represent that of an infi- 
nite box. Not only is sampling a key issue, but also the fact that 
simulations with periodic boundary conditions have no fluctu- 
ations on scales larger than the box size. To minimize and test 
for these effects we were conservative in our choices of box 
size and the mass range probed in each individual box. We 
also used nested-volume simulations to directly test for finite- 
volume effects. Because we used multiple boxes and aver- 
aged mass function results over the box ensemble, extended 
Press-Schechter theory can be used to correct for residual finit e 
volume-effects ("Mo & White1ll996l: iBarkana & Loebll2004 : 
this approach is differe nt from the individual box corrections 
applied by Reed et alT~l (|2007[) . Details are given in 35.31 

The paper is organized as follows. In ^ we give a brief 
overview of the mass function and popular fitting formulae, dis- 
cussing as well previous numerical work on the halo mass func- 
tion at high redshifts. In Owe give a short description of the 
A^-body code MC^ (Mesh-based Cosmology Code) and a sum- 
mary of the performed simulations. In ^we derive and discuss 
some simple criteria for the starting redshift and consider sys- 
tematic errors related to the numerical evolution such as mass 
and force resolution and time stepping. These considerations 
in turn specify the input parameters for the simulations in order 
to span the desired mass and redshift range for our investiga- 
tion. In Owe present results for the mass function at different 
redshifts as well as the halo growth function. Here we also dis- 



cuss the importance of post-processing corrections such as FOF 
particle sampling compensation and finite-volume effects. We 
discuss our results and conclude in O 

2. DEFINITIONS AND PREVIOUS WORK 

The mass function describes the number density of halos of 
a given mass. In order to determine the mass function in sim- 
ulations one has to first identify the halos and then define their 
mass. No precise theoretical basis exists for these operations. 
Nevertheless, depending on the situation at hand, the observa- 
tional and numerical communities have adopted a few "stan- 
dard" ways of defining halos and their associated masses. For 
a rec ent review o f these issues with regard to observations, see, 
e.g., 'Voit' ('2005), but for a more theoretically oriented review, 
see, e.g.. White (2001). 

2.1. Halo Mass 

There are basically two ways to find halos in a simulation. 
One, the overdensity method, is based on identifying overdense 
regions above a certain threshold. The threshold can be set with 
respect to the critical density pc = 3H^/8ttG (or the background 
density pb = ^mPc, where fl^ is the matter density of the Uni- 
verse including dark matter and baryons). The mass Ma of 
a halo identified this way is defined as the mass enclosed in 
a sphere of radius ta whose mean density is Apc- Common 
values for A range f rom 100 to 500 (or even higher). As ex- 
plained in I Voit I (2005), cluster observers prefer higher values 
for A. Properties of clusters are easier to observe in higher den- 
sity regions and these regions are more relaxed than the outer 
parts which are subject to the effects of inflow and incomplete 
mixing. The disadvantage of defining a halo in this manner is 
that sphericity of halos is implied, an assumption which may 
be easily violated, e.g., in the case of halos that formed in a re- 
cent merger event or halos at high redshifts. At higher redshifts, 
the nonlinear mass scale decreases rapidly, and the ratio of 
the considered halo mass Mhaio to can become large. This 
translates into producing large-scale structures roughly analo- 
gous to supercluster structures today. While these structures are 
gravitationally bound, they are often not virialized, nor spher- 
ical. Even the much smaller structures (which are considered 
in this paper) are not virialized at high redshifts, and therefore, 
assumptions about sphericity are most likely violated. Hence 
the spherical overdensity method does not suggest itself as an 
obvious way to identify halos at high redshift. 

The other method, the FOF algorithm, is based on finding 
neighbors of particles and neighbors of neighbors as defined by 
a given separation distance (see, e.g., Einasto et al. 1984; Davis 
et al. 1985). The FOF algorithm leads to halos with arbitrary 
shapes since no prior symmetry assumptions have been made. 
The halo mass is defined simply as the sum of particles which 
are members of the halo. While this definition is easy to apply 
to simulations, the connection to observations is difficult to es- 
tablish directly. (For an investigation of connections between 
different definitions of halos masses and approximate conver- 
sions between them, see White 2001). 

It is important to keep in mind that the definition of a halo is 
essentially the adoption of some sort of convention for the halo 
boundary. In reality, a sharp distinction between the particles 
in a halo and particles in the simulation "field" does not exist. 
Jenkins showed that the choice of a FOF finder with a linking 
length b = 0.2 to define halo masses provides the best fit for a 
universal form of the mass function. This choice has since been 
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adopted by many numerical practitioners as a standard conven- 
tion. A useful d i scussi on of the various halo definitions can be 
found in lWhitel(l200l . 

In this paper we use the FOF algorithm to identify halos and 
their masses. It was recently pointed out by Warren et al. (2006, 
hereafter Warren) that FOF masses suffer from a systematic 
problem when halos are sampled by relatively small numbers 
of particles. Although halos can be robustly identified with as 
few as 20 particles, if a given halo has too few particles, its 
FOF mass turns out to be systematically too high. We describe 
how we compensate for this effect in ^5.21 In the current paper, 
all results for the mass function are displayed at a fixed FOF 
linking length of b = 0.2, using the Warren correction. 

2.2. Defining the Mass Function 

The exact definition of the mass function, e.g., integrated ver- 
sus differential form or count versus number density, varies 
widely in the literature. To characterize different fits, Jenk- 
ins introduced the scaled differential mass function f{a,z) as 
a fraction of the total mass per Inci"' that belongs to halos: 

_ dp/pb M dn{M,z) 



fi'J,z) : 



(1) 



c/lncr-i pb(z)dln[a-^(M,z)]' 
Here n(M,z) is the number density of halos with mass M, pbiz) 
is the background density at redshift z, and a(M,z) is the vari- 
ance of the linear density field. As pointed out by Jenkins, this 
definition of the mass function has the advantage that to a good 
accuracy it does not explicitly depend on redshift, power spec- 
trum, or cosmology; all of these are encapsulated in a{M,z)- 
For the most part, we will display the mass function 

dn 

FiM,z)=-—— (2) 
fllogM 

as a function of logM itself. [In ^ we include results for 
f(<J,z).] 

To compute <j(M,z), the power spectrum P(k) is smoothed 
with a spherical top-hat filter function of radius R, which on 
average encloses a mass M(R = [3M/4Trpbiz)]^/^y. 

a^(M,z)=--^ / k^P(k)W\k,M)dk, (3) 



where W(k,M) is the top-hat filter: 



W(r) = 



W(k): 
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0,r>R 

[sm(kR)-kRcos(kR)]. 



(4) 



(5) 



(kR)' 

The redshift dependence enters only through the growth factor 
d{z), normalized so that d(Q) = 1 : 

a(M,z) = cr(M,0)d{z). (6) 
In the approximation of negligible difference in the CDM and 
baryon peculiar velocities, the growth function in a ACDM uni- 
verse is given by (Peebles 1980) 

D+(fl) 



d{a) -- 



(7) 



D+(a= 1)' 

where we consider d as a function of the cosmological scale 
factor a = 1/(1 +z), and 

5f2m Hia) f da' 



D^(a) = ■ 



2 Ho Jo [a'H(a')/Ho]^ 

I 1/2 



(8) 



with H(a)/Ho = [flm/a^ + Ci-- ^m)] ■ In particular, forz > 1, 
when matter dominates the cosmological constant, D'''(a) ~ a. 



Even in linear theory, equation ^ is only an approxima- 
tion because baryons began their gravitational collapse with 
velocities different from those of CDM particles. Until re- 
combination at z ^ 1100, well into the matter era with non- 
negligible growth of CDM inhomogeneities, the baryons were 
held against collapse by the p ressure of the CMB photons (see, 
e.g. IHu & Sugivama I (Il996h ). While thereafter the relative 
baryon-CDM velocity decayed as 1 /a, the residual velocity dif- 
ference was sufficient to affect the growth functi on d{z) at z = 50 
by niore than 1% and at z= 1 by about 0.2% dYoshida et al. I 
I200I lNaor& Barkana II200I . 

2.3. Fitting Functions 

Over the last three decades several different fitting forms for 
the mass function have been suggested. The mass function is 
not only a sensitive measure of cosmological parameters by it- 
self but also a key ingredient in analytic and semianalytic mod- 
eling of the dark matter distribution, as well as of several as- 
pects of the formation, evolution, and distribution of galaxies. 
Therefore, if a reliable and accurate fit for the mass function 
applicable to a wide range of cosmologies and redshifts were to 
exist, it would be of obvious utility. In this section we briefly 
review the common fitting functions and compare them at dif- 
ferent redshifts. 

The first analytic model for the mass function was developed 
by PS. Their theory accounts for a spherical overdense region 
in an otherwise smooth background density field, which then 
evolves as a Friedmann universe with a positive curvature. Ini- 
tially, the overdensity expands, but at a slower rate than the 
background universe (thus enhancing the density contrast), un- 
til it reaches the 'turnaround' density, after which collapse be- 
gins. Although from a purely gravitational standpoint this col- 
lapse ends with a singularity, it is assumed that in reality - due 
to the spherical symmetry not being exact - the overdense re- 
gion will virialize. For an Einstein-de Sitter universe, the den- 
sity of such an overdense region at the virialization redshift is 
z w 180pc(z)- At this point, the density contrast from the linear 
theory of perturbation growth [S(x,z) = d(z)S(x,Q)] would be 
^ciz) ~ 1.686 in an Einstein-de Sitter cosmology. For < 1, 
the value of the threshold parameter 4 can vary (see Lacey & 
Cole 1993), but the dependence on cosmology has little quan- 
titative significance (see, e.g., Jenkins). Thus, throughout this 
paper we adopt 6^ = 1 .686. 

Following the above reasoning and with the assumption that 
the initial density perturbations are described by a homoge- 
neous and isotropic Gaussian random field, the PS mass func- 
tion is specified by 



2 6c 

fps((j) = \ exp 



2^2 



(9) 



The PS approach assumes that all mass is inside halos, as en- 
forced by the constraint 



/+00 
/ps(a)Jlna-' = l. 
00 



(10) 



While as a first rough approximation the PS mass function 
agrees with simulations at z = reasonably well, it overpredicts 
the number of low-mass halos and underpredicts the number of 
massive halos at the current epoch. Furthermore, it is signifi- 
cantly in error at high redshifts (see, e.g., Springel et al. 2005; 
Heitmann et al. 2006a; ^SU l. 

After PS, several suggestions were made in order to improve 
the mass function fit. These suggestions were based on more 
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Table 1 
Mass Function Fits for f{a) 



Reference 



Fitting Function f{a) 



Mass Range 



Redshift range 



Reed et al. (2007) 



ST, Sheth & Tormen (2001) 0.3222 ^ 22 fexp 

Jenkins 0.315exp [-| Incr"' +0.61 p-^] 

Reed et al. (2003) /sT(cr) exp {-0.7/ [cr(cosh(2(T))5] } 

Warren 0.7234 (ct-' ^'^s + 0.2538) exp [-^^ 



1 + 



1 + 



x-^exp 



(^) +0.6Gi((T) + 0.4G2(a) 

0.03 ( <5, N "'^ 



ca&- 0.03 /'if 
20-2 (ncrf+3)- V cr / 



unspecified 
-1.2 <lncr-i < 1.05 
-1.7 < Incr-' < 0.9 
(lOi°-lO'^)/i-iM0 
-0.5 <lncr-' < 1.2 



unspecified 
z = 0-5 
z = 0-15 

z = 
z = 0-30 



Note. — Shown are examples of commonly used fitting functions. ST used a = 0.707 and p = 0.3, while I Sheth & Tormen 1 120021) suggest that a = 0.75 leads to a 
better fit. The Warren fit represents by far the largest uniform set of simulations based on multiple boxes with the same cosmology run with th e same code. We us e 
it as a reference standard throughout this paper. Reed et al. (2003) suggest an empirical adjustement of the ST fit, which is slightly modified in lReed et al. I i2007h . 
For the latter, Gi(cr) and G2(cr) are given by eqs. U6t and (TT), respectively, c = 1.08, ca = 0.764, and A = 0.3222. 



refined dynamical modeling, direct fitting to simulations, or a 
combination of the two. 

Using empirical arguments ST proposed an improved mass 
function fit of the form: 

^2 \ Pi 



fsT(a) = 0.3222^/ exp ( 



1 + 



(11) 



with a = 0.707 and p = 0.3. (Sheth & Tormen 2002 suggest 
a = 0.75 as an improved value.) Sheth et al. (2001) rederived 
this fit theoretically by extending the PS approach to an ellip- 
tical collapse model. In this model, the collapse of a region 
depends not only on its initial overdensity but also on the sur- 
rounding shear field. The dependence is chosen such tha t it re- 
covers the Zel'dovich approximation (IZel'dovich II 1 9701) in the 
linear regime. A halo is considered virialized when the third 
axis collapses (see also Lee & Shandarin (1998) for an earlier, 
different approach to the same idea). 

Jenkins combined high resolution simulations for four differ- 
ent CDM cosmologies (rCDM, SCDM, ACDM, and OCDM) 
spanning a mass range of over 3 orders of magnitude (~ (10'^- 
1O'^)/!~'M0), and including several redshifts between z = 5 
and 0. Independent of the underlying cosmology, the follow- 
ing fit provided a good representation of their numerical results 
(within ±20%): 

./je„ki„s('T) = 0.315exp(-|ln(T-'+0.61p'*) . (12) 

The above formula is very close to the Sheth-Tormen fit, lead- 
ing to some improvement at the high-mass end. The disadvan- 
tage is that it cannot be simply extrapolated beyond the range 
of the fit, since it was tuned to a specific set of simulations. 

By performing 16 nested-volume dark matter simulations. 
Warren was able to obtain significant halo statistics spanning 
a mass range of 5 orders of magnitude (~ (10'"- 1O'^)/i~'M0). 
Because this represents by far the largest uniform set of simulations- 
based on multiple boxes with the same cosmology run with the 
same code-we use it as a reference standard throughout this pa- 
per Using a functional form similar to ST, Warren determined 
the best mass function fit to be 

1.1982' 



/waiTen(cr) = 0.7234 (cr"' '^^^ + 0.2538) exp 



(13) 



For a quantitative comparison of the different fits at different 
redshifts, we show the ratio of the PS, Jenkins, and ST fits with 




7.5 



8 8.5 9 9.5 

log(M [MM) 



10 



Fig. 1 . — Ratio of the Jenkins, PS, and ST mass function fits with respect to 
the Warren fit for five different redshifts over a range of halo masses. Top to 
bottom: Redshifts z = 0, 5, 10, 15, and 20. Note that the ranges of the axes are 
different in the different panels. We do not show the Jenkins fit below masses 
of 10"/r'MQ at z = 0, since it is not valid for such low masses at that redshift. 



respect to the Warren fit in Figure[T] We do not show the Jenk- 
ins fit below 10" Ii'^Mq at z = since it diverges in this regime. 
The original ST fit, the Jenkins fit, and the Warren fit all give 
similar predictions. The discrepancy between PS and the other 
fits becomes more severe for higher masses at high redshifts. PS 
dramatically underpredicts halos in the high-mass range at high 
redshifts (assuming that the other fits lead to reasonable results 
in this regime). For low-mass halos the disagreement becomes 
less severe. For z = the Warren fit agrees, especially in the 
low-mass range below 10'^/!~'Mq, to better than 5% with the 
ST fit. At the high-mass end the difference increases up to 20%. 
The Jenkins fit leads to similar results over the considered mass 
range. At higher redshifts and intermediate-mass ranges around 
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10^ /z 'M0, the Warren and ST fit disagree by roughly a factor 
of 2. 

Several other groups have suggested modifications o f the ST 
fit. In^ we compare our results with two of them. Re ed et al. I 
(l2003h suggest an empirical adjustment to the ST fit by multi- 
plying it with an exponential function, leading to 

/Reed03(fT) = fsHa) cxp {-0.7/ [a(cosh(2a))'] } , (14) 

valid over the range -1.7 < In cr"' < 0.9. T his adjustment leads 
to a suppression of the ST fit at large . In lReed et al. I (l2007l) 
the adjustment to the ST fit is slightly modified again, leading 
to the following new fit: 



1000 



fReedOlicr) = Axl 

TT 




1+ ( ) +O.6G1+O.4G2 



0.03 



(neff+3)2 \ a 
ln(cr-i -0.4)2"! 



0.6 



2(0.6)2 
ln(cr-i -0.75)2 



(15) 

(16) 
(17) 



2(0.2)2 

with c = 1.08, ca = 0.764, and A = 0.3222. The a djustment has 
very s imilar effects to tha t of iReed et al. I ( l2003h . as we show 
in |5] iReed et al. I (l2007l) note that the (small) suppression of 
the mass function relative to ST as a function of redshift seen 
in simulations (see also Heitmann et al. 2006a) can be treated 
by adding an extra parameter, the power spectral slope at the 
scale of the halo radius, Weff (formally defined by equation (|42] | 
below). We return to this issue when we discuss our numerical 
results in §5. We summarize the described, most commonly 
used fitting functions in Table [1] 

Although fitting functions may be a useful way to approxi- 
mately encapsulate results from simulations, meaningful com- 
parisons to observations require overcoming many hurdles, e.g., 
an operational understanding of the definition of halo mass (see, 
e.g.. White 2001), how it relates to various observations, and 
error control in A^-body codes (see, e.g., O'Shea et al. 2005; 
Heitmann et al. 2005). In this paper, our focus is first on iden- 
tifying possible systematic problems in the A^-body simulations 
themselves and how they can be avoided and controlled. 

2.4. Halo Growth Function 

A useful way to study the statistical evolution of halo masses 
in simulations is to transform the mass functi on into the halo 
growth function, n(Mi , M2 , z) = F dlogM (iHeitmann et al. I 
l2006al) . which measures the mass-binned number density of ha- 
les as a function of redshift. The halo growth function, plotted 
versus redshift in Figure |2] shows at a glance how many halos 
in a particular mass bin and box volume are expected to exist 
at a certain redshift. This helps set the required mass and force 
resolution in a simulation which aims to capture halos at high 
redshifts. For a given simulation volume, the halo growth func- 
tion directly predicts the formation time of the first halos in a 
given mass range. 

In order to derive this quantity approximately, we first con- 
vert an accurate mass function fit (we use the Warren fit here) 
into a function of redshift z- It has been shown recently by us 
(Heitmann et al. 2006a) that mass function fits work reliably 
enough out to at least z = 20, and can therefore be used to esti- 
mate the halo growth function. Figure |2] shows the evolution of 
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Fig. 2. — Halo growth function based on the Warren mass function fit for 
different mass bins. The curves for the lower mass bins have a maximum at 
z > which reflects a crossover of the mass functions at different redshifts. 



eight different mass bins, covering the mass range investigated 
in this paper, as a function of redshift z- As expected from the 
paradigm of hierarchical structure formation in a ACDM cos- 
mology, small halos form much earlier than larger ones. An 
interesting feature in the lower mass bins is that they have a 
maximum at different redshifts. The number of the smallest ha- 
los grows until a redshift of z = 2 and then declines when halos 
start merging and forming much more massive halos. This fea- 
ture is reflected in a crossing of the mass functions at different 
redshifts for small halos. 

2.5. Mass Function at High Redshift: Previous Work 

Most of the effort to characterize, fit, and evaluate the mass 
function from simulations has been focused on or near the cur- 
rent cosmological epoch, z ^ 0. This is mainly for two reasons: 
( 1 ) so far most observational constraints have been derived from 
low -redshift objects (z < 1); (2) the accurate numerical evalua- 
tion of the mass function at high redshifts is a nontrivial task. 

The increasing reach of telescopes on the ground and in 
space, such as the upcoming James Webb Space Telescope, al- 
lows us to study the Universe at higher and higher redshifts. 
Recent discoveries include 970 galaxies at redshi fts be tween z = 

1.5 an d z = 5 from the VIMOS VLT Deep Survey dLe Fe vre et al~| 

'2005"), and the recent observation of a galaxy at z = 6.5 dMobasher et al. I 
2005). The epoch of reionization (FOR) is of central impor- 
tance to the formation of cosmic structure. Although our cur- 
rent observational knowledge of the FOR is rather limited, fu- 
ture 21 cm experiments have the potential for revolutionizing 
the field. Proposed low-frequency radio telescopes include LO- 
FAR (L ow Frequency Array) FL the Mileura Wide Field Array 
(MWA) dBowman et al. I l2006l2l and the next-generation SKA 
(Square Kilometer Array) 0. The observational progress is an 
important driver for high-redshift mass function studies. 



' See' http://www.lofar.org | 

^See http://haystack.mit.edu/arrays /MWA7] 

^ See http://www. skatelescope.org , 
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Table 2 



Summary of the Performed Runs 


Mesh 


DOa olZc 

(/i-'Mpc) 


ivcSUlULlUIl 

(/i-'kpc) 


Zin 




r oTLlClc iVldSS 


olllallcSL ridlO 


No. of Realizations 


1024^ 


256 


250 


100 





8.35 X 10'" 


3.34 X 10'^ 


5 


1024^ 


128 


125 


200 





1.04 X 10"' 


4.18 X 10" 


5 


1024^ 


64 


62.5 


200 





1.31 X 10*^ 


5.22 X 10'° 


5 


1024^ 


32 


31.25 


150 


5 


1.63 X 10*^ 


6.52 X 10'' 


5 


1024^ 


16 


15.63 


200 


5 


2.04 X 10^ 


8.16 X 10^ 


5 


1024^ 


8 


7.81 


250 


10 


2.55 X lO*" 


1.02 X 10^ 


20 


1024^ 


4 


3.91 


500 


10 


3.19 X 10^ 


1.27 X 10^ 


15 



Note. — Mass and force resolutions of the different runs. The smallest halos we consider contain 40 particles. All simulations have 256^^ particles. 



Theoretical studies of the mass function at high redshifts are 
challenging due to the small masses of the halos at early times. 
In order to capture these small-mass halos, high mass and force 
resolution are both required. For the large simulation volumes 
typical in cosmological studies, this necessitates a very large 
number of particles, as well as very high force resolution. Such 
simulations are very costly, and only a very limited number can 
be performed, disallowing exploration of a wide range of possi- 
ble simulation parameters. Alternatively, many smaller volume 
simulation boxes, each with moderate particle loading, can be 
employed. This leads automatically to high force and mass res- 
olution in grid codes (such as particle-mesh [PM]) and also re- 
duces the costs for achieving sufficient resolution for particle 
codes (such as tree codes) or hybrid codes (such as TreePM). 
The disadvantages of this strategy are the limited statistics in 
individual realizations (because fewer halos form in a smaller 
box) and the unreliability of simulations below an intermediate 
redshift at which the largest mode in the box is still (accurately) 
linear. In addition, results from small boxes may be biased, 
since they only focus on a small region and volume. Therefore, 
one must show that the simulations are free from finite-volume 
artifacts, e.g. missing tidal forces, and run a sufficient number 
of statistically independent simulations to reduce the sample 
variance. Both strategies, employing large volume or multi- 
ple small-volume simulations, have been followed in the past 
in order to obtain results at high redshifts. The different mass 
ranges investigated by different groups are shown in Figure [3] 
The fits are shown for redshifts z = 10 and 20. In the Appendix 
we provide a very detailed discussion on previous findings as 
organized by simulation volume. 

In summary, there is considerable variation in the high- 
redshift (z > 10) mass function as found by different groups, 
independent of box size and simulation algorithm. Broadly 
speaking, the results fall into two classes: either consistent with 
linear theory scaling of a universal form (Jenkins, Reed, ST, or 
Warren) at low redshift (Reed et al. 2003, 2007; Springel et 
al. 2005; Heitmann et al. 2006a; Maio et al. 2006; Zahn et 
al. 2007) or more consistent with the PS fit (Jang-Condell & 
Hernquist 2001; Yoshida et al. 2003a, 2003b, 2003c; Cen et al. 
2004; Iliev et al. 2006; Trac & Cen 2006). 

Our aim here is to determine the evolution of the mass func- 
tion accurately, at the few percent level, and at the same time 
characterize many of the numerical and physical factors that 
control the error in the mass f unction (details below) . We fol- 
low up on our previous work dHeitmann et al. Il2006a ) and ana- 



lyze a large suite of A^-body simulations with varying box sizes 
between 4 and 256/!~'Mpc, including many realizations of the 
small boxes, to study the mass function at redshifts up to z = 20 
and to cover a large mass range between 10^ and 10'^-^ /j-'Mq. 
With respect to our previous work, the number of small-box re- 
alizations has been increased to improve the statistics at high 
redshifts. Our results categorically rule out the PS fit as be- 
ing more accurate than any of the more modern forms at any 
redshift up to z = 20, the discrepancy increasing with redshift. 

3. THE CODE AND THE SIMULATIONS 

All simulations in this paper are carried out with the parallel 
PM code MC^. This code solves the Vlasov-Poisson equations 
for an expanding universe. It uses standard mass deposition and 
force interpolation methods allowing periodic or open bound- 
ary conditions with second-order (global) symplectic time step- 
ping and fast fourier transform based Poisson solves. Particles 
are deposited on the grid using the cloud-in-cell method. The 
overall computational scheme has proven to be accurate and 
efficient; relatively large time steps are possible with excep- 
tional energy conservation being achieved. MC^ has been ex- 
tensively tested against state-of-the-art cosmological simulation 
codes (Heitmann et al. 2005, 2007). 

We use the following cosmology for all simulations: 

r!,ot=1.0, r!cDM = 0.253, ribaiyon = 0.048, 



fT8 = 0.9, //o = 70kms"'Mpc"', n=l, 



(18) 



in concordance with cosmic microwave backgroun d and large 
scale structure observations jMacTavish et al. Il2006h (the third- 
year Wilkinson Microwav e Anisotropy P r obe ob servations sug- 
gest a lower value of erg; ISpergel et al. I ( 20071)). The transfer 
functions are generated with CMBFAST dSeliak & Zaidarriagal 
1996). We summarize the different runs, including their force 
and mass resolution, in Table |2] As mentioned earlier, we 
identify halos with a standard FOF halo finder with a linking 
length ofb = 0.2. Despite several shortcomings of the FOF halo 
finder, e.g., the tendency to link up two halos which are close 
to each other (see, e.g., Gelb & Bertschinger 1994, Summers et 
al. 1995) or statistical biases (Warren), the FOF algorithm itself 
is well defined and very fast. As discussed in ^2.11 we adopt the 
correction for sampling bias given by Warren when presenting 
our results. 

4. INITIAL CONDITIONS AND TIME EVOLUTION 

In a near-ideal simulation with very high mass and force res- 
olution, the first halos would form very early. By z = 50, a 
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Fig. 3. — Summary of recent work on the mass function at high redshift. 
The mass function fits are shown at z = 10 (top) and z = 20 (bottom) for the 
cosmology used throughout this paper (the other groups used slightly different 
parameters). At z = 10, Jang-Condell & Hemquist (2001) (gray shaded region) 
cover the very low mass range using a very small box, as do Cen et al. (2004) 
(green shaded region). The larger boxes of Reed et al. (2007) and Springel et 
al. (2005) (red shaded region) lead to results at higher halo masses. Note that 
in this regime the PS fit deviates substantially from the other fits, while at the 
very low mass end all fits tend to merge. Our suite of variable box sizes covers 
a mass range of 10^ to lO"^ /i^'Mq between z = and 20. a much larger range 
than previously covered by any group with a uniform set of simulations. At 
z = 20 Yoshida et al. (2003a, 2003b, 2003c, 2003d) cover the very low mass 
end of the mass function, while Zahn et al. (2007) investigate larger mass 
halos. Our simulations overlap with b oth of them at the edges. By combining 
a heterogeneous set of simulations. Re ed et al"! 420071) cover a wide range in 
mass and redshift. Figure quality reduced for the arXiv version of the paper. 



redshift commonly used to start cosmological simulations, a 
large number of small halos would already be present (see, e.g.. 
Reed et al. [2005] for a discussion of the first generation of 
star-forming halos). In a more reahstic situation, however, the 
initial conditions at z = 50 have of course no halos, the particles 



le+06 I 




Fig. 4. — Probability distribution of | V</)[ in units of the interparticle spacing 
Ap. All curves shown are drawn from 256'' particle simulations from an initial 
density grid of 256' zones. The physical box sizes are 126/r'Mpc (black 
line), 32/r'Mpc (red line), and 8/i"'Mpc (green line). As expected, (|V(/)|) 
increases with decreasing box size (which is equivalent to increasing force 
resolution). Therefore, Zm snd Zcmss are higher for the smaller boxes. 



having moved only the relatively small distance assigned by the 
initial Zel'dovich step. Only after the particles have traveled a 
sufficient distance and come close together can they interact lo- 
cally to form the first halos. In the following we estimate the 
redshift when the Zel'dovich grid distortion equals the interpar- 
ticle spacing, leading to the most conservative estimate for the 
redshift of possible first halo formation. From this estimate, we 
derive the necessary criterion for the starting redshift for a given 
box size and particle number. 

4.1. Initial Redshift 

In order to capture halos at high redshifts, we have found 
that it is very important to start the simulation sufficiently early. 
We consider two criteria for setting the starting redshift: (1) en- 
suring the linearity of all the modes in the box used to sample 
the initial matter power spectrum, and (2) restricting the initial 
particle move to prevent interparticle crossing and to keep the 
particle grid distortion relatively small. The first criterion is 
commonly used to identify the starting redshift in simulations. 
However, as shown below, it fails to provide sufficient accuracy 
of the mass functions, accuracy which can be obtained when 
a second (much more restrictive) control is applied. Further- 
more, it is important to allow a sufficient number of expansion 
factors between the starting redshift zin and the highest redshift 
of physical significance. This is needed to make sure that arti- 
facts from the Zel'dovich approximation are negligible and that 
the memory of the artificial particle distribution imposed at Zin 
(grid or glass) is lost by the time any halo physics is to be ex- 
tracted from the simulation results. 

Although not studied here, it is important to note that high- 
redshift starts do require the correct treatment of baryons as 
noted in ^2.21 In addition, redshift starts that are too high can 
lead to force errors for a variety of reasons, e.g., interpolation 
systematics, round-off, and correlated errors in tree codes. 
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Table 3 

Initial Redshift Estimates from the Linearity of 

A'(fcNy) 



Box Size 
(/j-'Mpc) 


(/iMpc"') 


r(z=0,^Ny) 




126 


6.3 


0.0002 


33 


32 


25 


1.7-10-5 


45 


16 


50 


4.8-10"^ 


50 


8 


100 


1.3-10"^ 


55 



Note. — The number of particles is 256^^, tlie same in all simulations. 



4.1.1. Initial Perturbation Amplitude 

The initial redshift in simulations is often determined from 
the requirement that all mode amplitudes in the box below 
the particle Nyquist wavenumber characterized by with 
kf^y = 27r/Ap, where Ap is the mean interparticle spacing, be 
sufficiently linear. The smaller the box size chosen (keeping the 
number of particles fixed), the larger the largest A;-value. There- 
fore, in order to ensure that the smallest initial mode in the box 
is well in the linear regime, the starting redshift must increase 
as the box size decreases. In the following we give an estimate 
based on this criterion for the initial redshift for different sim- 
ulation boxes. We (conservatively) require the dimensionless 
power spectrum A^ = k^P(k)/27T^ to be smaller than 0.01 at the 
initial redshift. The initial power spectrum is given by 

A(^Ny,Zi„)=^^^ 2,2(,,„+l)2 ' (19) 

where B is the normalization of the primordial power spectrum 
(see, e.g., Bunn & White [1997] for a fitting function for B in- 
cluding COBE results) and T{k) is the transfer function. We 
assume the spectral index to be n = 1, which is sufficient to ob- 
tain an estimate for the initial redshift. For a ACDM universe 
the normalization is roughly B 3.4 x lQ^{h~^Mpc)'^ . There- 
fore, Zin is simply determined by 

Zi„~4150 4yr(z = 0,^Ny)- (20) 
We present some estimates for different box sizes in Table [3] 
For the smaller boxes (< 8 /i-'Mpc), the estimates for the initial 
redshifts are at around Zin = 50. 

It is clear that this criterion simply sets a minimal require- 
ment for Zin and neglects the fact that the initial particle move 
should be small enough to maintain the dynamical accuracy of 
perturbation theory (linear or higher order) used to set the ini- 
tial conditions. Also, this criterion certainly does not tell us 
that if, e.g., Zin = 50, then we may already trust the mass func- 
ti on at, say, z = 30. An example of this is provided by the results 
of iReed et al.1 ( l2003l) . who find that their high-redshift results 
between z = 7 and 15 have not converge if they start their sim- 
ulations at Zin = 69. (A value of Zin =139 was claimed to be 
sufficient in their case.) 

We now consider another criterion - ostensibly similar in 
spirit - that particles should not move more than a certain frac- 
tion of the interparticle spacing in the initialization step. This 
second criterion demands much higher redshift starts. 

4.1.2. First Crossing Time 

In cosmological simulations, initial conditions are most of- 
ten generated using the Zel'dovich approximation (.Zel'dovichl 




1 10 100 

Box length [Mpc/h] 



400p-rrf 




1 10 100 

Box length [Mpc/h] 



Fig. 5. — Average redshift of first crossing (top) and highest redshift of 
first crossing (bottom) as a function of box size. The initial conditions (five 
different realizations) are shown for boxes between 1 and 512/i^'Mpc with 
128^ and 256' particles. For each initial condition, zj?,"', and z™^ are shown 
by the crosses. The solid lines show the average from the five realizations. 
As expected, scatter from the different realizations is larger for smaller boxes. 
These plots provide estimates of the required initial redshift for a simulation 
since | V(}!>|/Ap is z-independent in the Zel'dovich approximation (see text). 

Il970l) . Initially each particle is placed on a uniform grid or in a 
glass configuration and is then given a displacement determined 
by the relation 

x = q-c/(z)V0, (21) 

where q is the Lagrangian coordinate of each particle. The gra- 
dient of the potential ((> is independent of the redshift z. The 
Zel'dovich approximation holds in the mildly nonlinear regime, 
as long as particle trajectories do not cross each other (no caus- 
tics have formed). Studying the magnitude of | V0| allows us to 
estimate two important redshift values: first, the initial redshift 
Zin at which the particles should not have moved on average 
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more than a fraction of the interparticle spacing Ap = Lbox/wp, 
where Lbox is the physical box size and «p the number of parti- 
cles in the simulation; second, the redshift at which particles 
first move more than the interparticle spacing, Zcioss, i-c, at 
which they have traveled on average a distance greater than Ap. 

For a given realization of the power spectrum, the magnitude 
of |V0| depends on two parameters: the physical box size and 
the interparticle spacing. Together these parameters determine 
the range of scales under consideration. The smaller the box, 
the smaller the scales; therefore, |V(/)| increases and both Zin 
and Zcross increase. Increasing the resolution has the same ef- 
fect. In Figure |4] we show the probability distribution function 
for |V0| for three different box sizes, 8, 32, and 126/!~'Mpc, 
representing values studied by other groups, as well as in this 
paper. To make the comparison between the different box sizes 
more straightforward, we have scaled |V0| with respect to the 
interparticle spacing Ap. All curves are drawn from simula- 
tions with 256^ particles on a 256^ grid, in accordance with the 
set up of our initial conditions. The behavior of the probability 
function follows our expectations: the smaller the box, or the 
higher the force resolution, the larger the initial displacements 
of the particles on average. From the mean and maximum val- 
ues of such a distribution we can determine appropriate values 
for Zin and Zcross- For our estimates we assume d(z) — 1/(1 
which is valid for high redshifts. The maximum and rms initial 
displacements of the particles can then be easily calculated: 

,m.x max(|V0|/Ap) 



1+Zin 

rms(|V(/)|/Ap) 

1+Zin 



(22) 



(23) 



The very first "grid crossing" of a particle occurs when (Jj™" = 
1 ; on average the particles have moved more than one particle 
spacing when = 1 . This leads to the following estimates: 

zfi;:L-max(V(/./Ap)-l, (24) 
CoL-rms(V0/Ap)-l. (25) 

We show these two redshifts in Figure |5] for 10 different box 
sizes ranging from 1 to 5 12 /!"'Mpc and for 256-^ and 128^* parti- 
cles. The top panel shows the average redshift of the first cross- 
ing as a function of box size (which corresponds to the max- 
imum in Fig. m. The bottom panel shows the redshift where 
the first "grid crossing" occurs (corresponding to the right tail 
in Fig. m. To estimate the scatter in the results, we have gen- 
erated five different realizations for each box. As expected, the 
small boxes show much more scatter. The average redshift of 
the first crossing in the 1 /z"^Mpc box varies between z = 63 
and 83, while there is almost no scatter in the 512/i"'Mpc box. 
Since | V0|/Ap is independent of redshift in the Zel'dovich ap- 
proximation, a simple scaling determines the appropriate initial 
redshift from these plots. For example, if a particle should not 
have moved more than 0.3Ap on average at the initial redshift, 
the average redshift of first crossing has to be multiplied by a 
factor 1/0.3 = 3.3. For an 8/!~'Mpc box this leads to a mini- 
mum starting redshift of z = 230, while for a 126/;~'Mpc box 
this suggests a starting redshift of Zin = 50. The 128^ particle 
curve can be scaled to the 256^ particle curve by multiplying 
by a factor of 2. Curves for different particle loadings can be 
obtained similarly. 

4.2. Transients and Mixing 



The Zel'dovich approximation matches the exact density and 
velocity fields to linear order in Lagrangian perturbation the- 
ory. Therefore, there is in principle an error arising from the 
resulting discrepancy with the density and velocity fields given 
by the exact growing mode initialized in the far past. 

This error is linear in the number of expansion factors be- 
tween Zin and the redshift of interest Zph vs- It has been e xplored 
in the context of simulation error by Valag easl (12002') and by 
ICrocce et al. I ([2006) ■ Depending on the quantity being cal- 
culated, the number of expansion factors between Zin and Zphys 
required to limit the error to some given value may or may not 
be easy to estimate. For example, unlike quantities such as the 
skewness of the density field, there is no analytical result for 
how this error impacts the determination of the mass function. 
Neither does there exist any independent means of validating 
the result aside from convergence studies. Nevertheless, it is 
clear that to be conservative, one should aim for a factor of ~ 20 
in expansion factor in order to anticipate errors at the several 
percent level, a rule of thumb that has been followed by many 
A^-body practitioners (and often violated by others!). This rule 
of thumb gives redshift starts that are roughly in agreement with 
the estimates in the previous subsection. Convergence tests 
done for our simulations show that the suppression in the mass 
function is very small (less than 1 %) for simulations whose evo- 
lution covers a factor of 15 in the expansion factor and can be 
up to 20% for simulations that evolved by only 5 expansion 
factors. However, due to modest particle loads, we were unable 
to distinguish between the error induced by too few expansion 
factors and the breakdown of the Zel'dovich approximation. 

Another possible problem, independent of the accuracy of 
the Zel'dovich approximation, is the initial particle distribu- 
tion itself. Whether based on a grid or a glass, the small- 
distance (k > kj^y) mass distribution is clearly not sampled at 
all by the initial condition. Therefore, unlike the situation 
that would arise if a fully dynamically correct initial condition 
were given, some time must elapse before the correct small- 
separation statistics can be established in the simulation. Thus, 
all other things being equal, for the correct mass function to 
exist in the box, one must run the simulation forward by an 
amount sufficiently greater than the time taken to establish 
the correct small-scale power on first-halo scales while eras- 
ing memory on these scales of the initial conditions. If this is 
not done, structure formation will be suppressed, leading to a 
lowering of the halo mass function. 

Because there is no fully satisfactory way to calculate Zin in 
order to compute the mass function at a given accuracy, we sub- 
jected every simulation box to convergence tests in the mass 
function while varying Zin- The results shown in this paper are 
all converged to the sub-percent level in the mass function. We 
give an example of one such convergence test below. 

4.2.1. Initial Redshift Convergence Study 

As mentioned above, we have tested and validated our esti- 
mates for the initial redshift for all the boxes used in the sim- 
ulation suite via convergence studies. Here, we show results 
for an 8 /!"'Mpc box with initial redshifts Zin = 50, 150, and 250 
in Figure |6l where the mass functions at z = 10 are displayed. 
For the lowest initial redshift, Zin = 50, the average initial par- 
ticle movement is 1.87Ap, while some particles travel as much 
as 5.03Ap. This clearly violates the requirement that the initial 
particle grid distortion be kept sufficiently below 1 grid cell. 
The starting redshift Zin =150 leads to an average displacement 




Fig. 6. — Dependence of the mass function on the initial redshift. The results are at z = 10 from thi'ee 8/1 'Mpc box simulations with zin = 50 (left), zin = 150 
(middle), and at ziR = 250 (right). The mass function in the left panel is systematically lower than the other two by roughly 15%. Poisson error bars are shown. 



of 0.63Ap and a maximum displacement of 1.71 Ap, and there- 
fore just barely fulfills the requirements. For Zm = 250 we find 
an average displacement in this particular realization of 0.37Ap 
and a maximum displacement of l.OOAp. 

The bottom plot in each of the three panels of Figure|6]shows 
the ratio of the mass functions with respect to the Warren fit. 
In the middle and right panels the ratio for the largest halo is 
outside the displayed range. The mass function from the simu- 
lation started at Zin = 50 (left panel) is noticeably lower, 15%, 
than for the other two simulations. The mass functions from the 
two higher redshift starts are in good agreement, showing that 
the choice for average grid distortion of approximately 0.3Ap 
is conservative, and that one can safely use (0.5-0.6)Ap. The 
general conclusion illustrated by Figure|6]is that if a simulation 
is started too late, halos are found to be missing over the entire 
mass range. With the late start, there is less time to form bound 
objects. Also, some particles that are still streaming towards a 
halo do not have enough time to join it. Both of these artifacts 
lead to an overall downshift of the mass function. 

To summarize, requiring a limit on initial displacements sets 
the starting redshift much higher than simply demanding that 
all modes in the box stay linear Indeed, the commonly used 
latter criterion (with 5™'* ~ 0.1) is not adequate for computing 
the halo mass function at high redshifts. One must verify that 
the chosen Zin sets an early enough start as shown here. We 
comment on previous results from other groups with respect to 
this finding below in ^ 

4.3. Force and Mass Resolution 

We now take up an investigation of the mass and force reso- 
lution requirements. The first useful piece of information is the 
size of the simulation box: from Figure |2] we can easily trans- 
late the number density into when the first halo is expected to 
appear in a box of volume V. For example, a horizontal line at 
n = 10"^ would tell us at what redshift we would expect on aver- 
age to find 1 halo of a certain mass in a (100/i~'Mpc)-^ box. The 
first halo of mass 10" - W^Iz-^Mq will appear at z ~ 15.5, 
and the first cluster-like object of mass 10'"*- lO'^/j'^M© at 
z ~ 2. Of course, these statements only hold if the mass and 
force resolution are sufficient to resolve these halos. The mass 
of a particle in a simulation, and hence the halo mass, is deter- 
mined by three parameters: the matter content of the Universe 



including baryons and dark matter, the physical box size 
, and the number of simulation particles n^: 



«VrtMe = 2.775 X 10"f^„ 



nph 



Mpc 



(26) 



The required force resolution to resolve the chosen smallest 
halos can be estimated very simply. Suppose we aim to resolve 
a virialized halo with comoving radius ta at a given redshift z, 
where A is the overdensity parameter with respect to the critical 
density pc- The comoving radius ta is given by 



rA=9.51 X 10" 





"1 


^1 Mac \ 






vA/i-iMoj 



1/3 



/r'Mpc, (27) 



where Q,{z) = \ +z)^/[rim(l +z)'' + i^A] and the halo mass 
Mac = nipaanh, where is the number of particles in the halo. 
We measure the force resolution in terms of 



5 ~ 



(28) 



In the case of a grid code, Wg is literally the number of grid 
points per linear dimension; for any other code, n„ stands for the 
number of "effective softening lengths" per linear dimension. 
To resolve halos of mass Mac, a minimal requirement is that 
the code resolution be smaller than the radius of the halo we 
wish to resolve: 

6f<rA. (29) 

Note that this minimal resolution requirement is aimed only 
at capturing halos of a certain mass, not at resolving their in- 
terior profile. Next, inserting the expression for the particle 
mass (eq. 1261 ) and the comoving radius (eq. ll27l ) into the 
requirement (eq ||29l ) and employing the relation between the 
interparticle spacing Ap and the box size Ap = Lbox/wp, the res- 
olution requirement reads 



Sf 

IT- < 0.62 
A„ 



n/,ri(z) 



A 



1/3 



(30) 



We now illustrate the use of this simple relation with an ex- 
ample. Let A = 200 and consider a ACDM cosmology with 
ftm = 0.3. Then for PM codes for which (5f/Ap = Wp/wg, we 
have the following conclusions. If the number of mesh points 
is the same as the number of particles (wp = Wg), halos with less 
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log(M [Wo/h]) log(M [ife/h]) loq(M [«©/h]) 



Fig. 7. — Convergence of the mass function as a function of force resolution. All results are shown at z = 0, for 256^ particles and a 126/i"'Mpc box with Poisson 
eiTor bars. The resolution varies between 256'' (left), 512' (middle), and 1024'' grid points (right). The vertical line denotes the predicted theoretical resolution 
limit: halos on the right of the line should not be lost. The resolution limit is 2500 particles per halo for the 256' giid, 300 particles per halo for the 512' grid, and 
40 particles per halo for the 1024' grid. 



than 2500 particles cannot be accurately resolved. If the num- 
ber of mesh points is increased to 8 times the particle num- 
ber (np = 1 /2ng), commonly used for cosmological simulations 
with PM codes, the smallest halo reliably resolved has roughly 
300 particles, and if the resolution is increased to a ratio of 1 
particle per 64 grid cells, which we use in the main PM sim- 
ulations in this paper, halos with roughly 40 particles can be 
resolved. It has been shown in Heitmann et al. (2005) that this 
ratio (1:64) does not cause collisional effects and that it leads to 
consistent results in comparison to high-resolution codes. Note 
that increasing the resolution beyond this point will not help, 
since it is unreliable to sample halos with too few particles. 
Note also that a similar conclusion holds for any simulation 
algorithm and not just for PM codes. 

In Figure [T] we show results from a resolution convergence 
test at z = 0. We run 256^ particles in a 126/i~'Mpc box with 
three different resolutions: 0.5, 0.25, and 0.125 /i~'Mpc. The 
vertical line in each figure shows the mass below which the 
resolution is insufficient to capture all halos following condi- 
tion ( [30l ). In all three cases, the agreement with the theoretical 
prediction is excellent. 

4.4. Time Stepping 

Next, we consider the question of time-step size and estimate 
the minimal number of time steps required to resolve the halos 
of interest. We begin with a rough estimate of the characteristic 
particle velocities in halos. For massive halos, the halo mass 
M20() and its velocity dispersion are connected by the approxi- 
mate relation (Evrard 2004): 



M200 ' 



(31) 



H/Ho Vl080km/s, 

A more accurate expression can be found in lEvrard et al. I (l2007l) . 
but the above is more than sufficient for our purposes. At high 
redshift, JIa can be neglected, and we can express the velocity 
dispersion as a function of redshift: 



M200 



1/3 



kms 



(32) 



In a time St, the characteristic scale length 61 is given by SI : 
(TySt or 



St' 



61 



100J//km / M200 



-1/3 



(33) 



The scale factor a is a convenient time variable for codes work- 
ing in comoving units, such as ours. Expressed in terms of the 
scale factor, equation (33[ reads: 



da ~ lO"* 



61 



M200 



-1/3 



(34) 



/!"'Mpc \h-^MQ^ 

We are interested in the situation where 61 is actually the force 
resolution, 6f. In a single time step, the distance moved should 
be small compared to 6f; i.e., the actual time step should be 
smaller than 6a estimated from the above equation when 61 
is replaced on the right-hand side with 6f. Let us consider 
a concrete example for the case of a PM code where 6{ = 
Lbox/no as explained earlier. For a "medium" box size of Lbox = 
256/!"'Mpc and a grid size of «g = 1024, = 0.25/i"'Mpc. For 
a given box, the highest mass halos present have the largest 
<jy and give the tightest constraints on the time step. For the 
chosen box size, a good candidate halo mass scale is M200 ^ 
10'^ H'^Mq (this could easily be less, but it does not change the 
result much). In this case, 

(5fl~ 0.025. (35) 

If, for illustration, we start a simulation at z = 50 and evolve 
it down to z = 0, this translates to roughly 40 time steps. We 
stress that this estimate is aimed only at avoiding disruption of 
the halos themselves, and is certainly not sufficient to resolve 
the inner structure of the halo. 

In Figure [8] we show two tests of the time step criterion. The 
top panel shows the result from a 32/!"'Mpc box at redshift 
z = 5. The simulation starts at Zin =150 and is evolved with 50, 
125, and 250 time steps down to z = 5. Following the argument 
above for this box size, one would expect all three choices to 
be acceptable, and the excellent agreement across these runs 
testifies that this is indeed the case. We also carried out a run 
with only five time steps, which yields a clearly lower (~ 20%) 
mass function than the others, but not as much as one would 
probably expect from such an imprecise simulation. 
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The bottom panel shows the results from a 126/i~^Mpc box 
at z = 0. This simulation was started at zin = 50 and run to 
Z = with 5, 8, 100, and 300 time steps. Again, as we would 
predict, the agreement is very good for the last two simulations, 
and the convergence is very fast, confirming our estimate that 
only 0(10) time steps is enough to get the correct halo mass 
function. Overall, the halo mass function appears to be a very 
robust measure, not very sensitive to the number of time steps. 
Nevertheless, we used a conservatively large number of time 
steps, e.g., 500 for the simulations stopping at z = and 300 for 
those stopping at z = 10. 




1 1 

log(M [M@/h]) 




log(M [M@/t)]) 



Fig. 8.— Top: One of the 32/!"'Mpc box realizations ran with 250, 125, 50 
and 5 time steps between zi„ =150 and zflnai = 5. The mass function is shown 
at the final redshift z = 5. Data points for all runs except the one with five time 
steps are so close that they are difficult to distinguish. Bottom: A 126/r']yipc 
box with 300, 100, 8, and 5 time steps between Zm — 50 and Zfinal = 0. The 
agreement for the very large halos for 100 and 300 time steps is essentially 
perfect. Poisson en'or bars are shown. 



In the previous subsections we have discussed and tested dif- 
ferent error control criteria for obtaining the correct simulated 
mass function at all redshifts. These criteria are (1) a suffi- 
ciently early starting redshift to guarantee the accuracy of the 
Zel'dovich approximation at that redshift and provide enough 
time for the halos to form; (2) sufficient force and mass resolu- 
tion to resolve the halos of interest at any given redshift; and (3) 
sufficient numbers of time steps. Violating any of these criteria 
always leads to a suppression of the mass function. Most sig- 
nificantly, our tests show that a late start (i.e., starting redshift 
too low) leads to a suppression over the entire mass range under 
consideration, and is a likely explanation of the low mass func- 
tion results in the literature. As intuitively expected, insufficient 
force resolution leads to a suppression of the mass function at 
the low-mass end, while errors associated with time stepping 
are clearly subdominant and should not be an issue in the vast 
majority of simulations. 

5. RESULTS AND INTERPRETATION 

In this section we present the results from our simulation 
suite. We describe how the data are obtained as well as the 
post-processing corrections applied. The latter include com- 
pensation for FOF halo mass bias induced by finite (particle 
number) sampling, and the (small) systematic suppression of 
the mass function induced by the finite volume of the simula- 
tion boxes. 

5.1. Binning of Simulation Data 

Before venturing into the simulation results, we first describe 
how they were obtained and reported from individual simula- 
tions. We used narrow mass bins while conservatively keep- 
ing the statistical shot noise of the binned points no worse than 
some given value. Bin widths AlogM were chosen such that 
the bins contain an equal number of halos A'h. The worst- 
case situation occurs at z = 20 for the 8 /i"'Mpc box, which has 
A^h = 80; the 4/!"'Mpc box at the same redshift has A^h = 400. At 
z= 15 we have A^h = 150, 1600, and 3000 for box sizes 16, 8, and 
4/i"'Mpc, respectively. At z = 10 the smallest value A'h = 450 
is for the 32/!"'Mpc box, while at z = 5 and we essentially 
always have A'h > 10000. 

With a mass function decreasing monotonically with M, this 
binning strategy results in bin widths increasing monotonically 
with M. The increasing bin size may cause a systematic de- 
viation - growing towards larger masses - from an underlying 
"true" continuous mass function. The data points for the binned 
mass function give the average number of halos per volume in 
a bin, 

F^A'hAVAlogM), (36) 

plotted versus an average halo mass, averaged by the number of 
halos in the bin: 

M=^M/A'h. (37) 

bin 

Assuming that the true mass function dn/d logM has some ana- 
lytic form F{M), a systematic deviation due to the binning pre- 
scription 

F-Fm ^33^ 



Ebir 



F{M) 



can be evaluated by computing F and M as 



F = 



AlogM' 



M = 



(39) 
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where dn = F{M)d\ogM and the integrations are over a mass 
range [M,M + AM]. For the leading-order term of the Taylor 
expansion of ebi„(Z^M), we find 

F"-2{F'f/F 



Cbin 



2AF 



(40) 



where the primes denote d/dM. A characteristic magnitude of 
this ebin for a general F{M) is (AM/M)^/24. However, in our 
case, where the relevant scales k':^ k^,^^ 0.01/;Mpc"', ebin has 
a much stronger suppression, as explained below. 

We know that the mass function is close to the universal 
form, 

Pb „, , d\na~^ 



F{M)=—f(a)- 



(41) 



(see, eq. |[T]). Note that for k ^ k^.^^, (T~'(M) is a slowly varying 
function, i.e., 

t/log(T~' _ neff+3 



dlogM 



(42) 



is much smaller than unity, and the derivative (ilog(T"'/c/logM 
also changes slowly with M. Then, despite the steepness 
of F{a) at small a, the factor /((T)iiln(T"'/'^logM in equa- 
tion (|4TJ depends weakly on M. Therefore, the mass function 
F{M) is close to being inversely proportional to M. In the limit 
of exact inverse proportionality, F cx M"', equation (|40] | tells 
us that ebin 0. This effective cancellation of the two terms 
on the right-hand side of equation (|40] | makes the binning error 
negligible to the accuracy of our F{M) reconstruction when- 
ever a bin width AlogM does not exceed 0.5. To confirm the 
absence of any systematic offsets due to the binning, we binned 
the data into logM intervals 5 times narrower and wider, with 
no apparent change in the inferred F(M) dependence. 

We remark that the situation could be quite different with 
another binning choice. For example, if the binned masses M 
were chosen at the centers of the corresponding logM intervals, 
logM = [logM+log(M+ AM)]/2, the systematic binning devi- 
ation 

' -(AM)2 (43) 



(center) 
^bin ' 



24F 



would have no special cancellation for the studied type of mass 
function. A corresponding binning error would be about 2 or- 
ders of magnitude larger than that of equations (|36] | and dJTl i. 

The statistical error bars used are Poisson errors, following 
the improved definition of Heinrich (2003): 



(44) 



At large values of Nh, these error bars asymptote to the familiar 
form v^h- At smaller values of A'h - which are of minor con- 
cern here - equation |44] has several advantages over the stan- 
dard Poisson error definition, some being (1) it is nonzero for 
A'h = 0; (2) the lower edge of the error bar does not go all the 
way to zero when A'h = 1 ; (3) the asymmetry of the error bars 
reflects the asymmetry of the Poisson distribution. 

Finally, as noted earlier and discussed in the next section, 
all the results shown in the following include a correction for 
the sampling bias of FOF halos according to equation (05]). 
This mass correction brings down the low-mass end of the mass 
function. 
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Fig. 9. — FOF mass con'ection for halos in 4 (dark blue), 8 (black), 16 
(light blue), and 32 (yellow) ft^'Mpc boxes. To show the effect clearly, we 
plot the ratio of our data to the Wan'en fit. Crosses show the unconnected mas s 
function and squares the mass function after connection, following eq. (45). 
Note the smooth behavior of the conected mass function as opposed to the 
mass-function jumps across box sizes for the unconected data. 



5.2. FOF Mass Correction 

The mass of a halo as determined by the FOF algorithm dis- 
plays a systematic bias with the number of particles used to 
sample the halo. Too few particles lead to an increase in the es- 
timated halo mass. By systematically subsampling a large halo 
population from N-body simulations (at z = 0), Warren deter- 
mined an empirical correction for this undersampling bias. For 
a halo with nh particles, his correction factor for the FOF mass 
is given by 

«r=«h(i-«h"-'). (45) 



We have carried out an independent exercise to check the sys- 
tematic bias of the FOF halo mass as a function of particle num- 
ber based on Monte Carlo sampling of an NFW halo mass pro- 
file with varying concentration and particle number, as well as 
by direct checks against simulations (e.g.. Fig. |9]l; our results 
are broadly consistent with equation ( l45T l. Details will be pre- 
sented elsewhere (Z. Lukic et al., in preparation). In this asso- 
ciated work we also address how overdensity masses connect to 
FOF masses, how this relation depends on the different linking 
length used for the FOF finder, and the properties of the halo 
itself, such as the concentration. 

The effect of the FOF sampling correction can be quickly 
gauged by considering a few examples: for a halo with 50 par- 
ticles, the mass reduction is almost 10%, for a halo with 500 
particles, it is ^ 2.4%, and for a well-sampled halo with 5000 
particles, it is only 0.6%. As a cautionary remark, this correc- 
tion formula does not represent a general recipe but can depend 
on variables such as the halo concentration. Since the condi- 
tions under which different simulations are carried out can dif- 
fer widely, corrections of this type should be checked for appli- 
cability on a case-by-case basis. Note also that the correction 
for the mass function itself depends on how halos move across 
mass bins once the FOF correction is taken into account. 

The choice of the mass function range in a given simula- 
tion box always involves a compromise: too wide a dynamic 
range leads to poor statistics at the high-mass end and possible 
volume-dependent systematic errors, and too narrow a range 
leads to possible undersampling biases. Our choice here re- 
flects the desire to keep good statistical control over each mass 
bin at the expense of wide mass coverage, compensating for 
this by using multiple box sizes. Therefore, in our case it is 
important to demonstrate control over the FOF mass bias. An 
example of this is shown in Figure |9] where results from four 
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box sizes demonstrate the successful application of the Warren 
correction to simulation results at z = 10. 

5.3. Simulation Mass and Growth Function 

The complete set of simulations, summarized in Table |2] al- 
lows us to study the mass function spanning the redshift range 
from z = 20 to 0. The mass range covers dwarf to massive 
galaxy halos at z = (cluster scales are best cov ered by much 
bigger boxes as in Warren and lReed et al. 1120071) . and at higher 
redshifts goes down to 10^/z~'Mq, the mass sc ale above which 
gas in halos can cool via atomic line cooling (iTegmark et al. I 
Il997h . 

5.4. Time Evolution of the Mass Function 

Halo mass functions from the multiple-box simulations are 
shown in Figure [10] with results being reported at five differ- 
ent redshifts with no volume corrections applied. The combi- 
nation of box sizes is necessary because larger boxes do not 
have the mass resolution to resolve very small halos at early 
redshifts, while smaller boxes cannot be run to low redshifts. 
The bottom plot of each panel shows the ratio of the numeri- 
cally obtained mass function, and various other fits, to the War- 
ren fit as scaled by linear theory (for volume-corrected results, 
see Fig.[T2]i. Displaying the ratio has the advantage over show- 
ing relative residuals that large discrepancies (more than 100%) 
appear more clearly. For all redshifts, the agreement with the 
Warren fit is at the 20% level. The ST fit matches the simula- 
tions for small masses very well but overpredicts the number 
of halos at large masses. This overprediction becomes worse at 
higher redshifts. For example, at z = 15 ST overpredicts halos 
of IO'^/z-'Mq by a factor of 2. Reed et al. (2003) found a sim- 
ilar result: the ST fit at z = 15 for halos with mass larger than 
lO'^r'M© disa grees with their simulation by 50%. Agree- 
ment with the Reed et al. (2003, 2007) fits is also good, within 
the 10% level. (For a further discussion focused around the 
question of universality, see Section 5.7.) The PS fit in general 
is not satisfactory over a larger mass range at any redshift. It 
crosses the other fits at different redshifts for different masses. 
Away from this crossing region, however, the disagreement can 
be as large as an order of magnitude, e.g. for z = 20 over the 
entire mass range we consider here. 

5.5. Halo Growth Function 

As discussed in 32.41 the halo growth function (the number 
density of halos in mass bins as a function of redshift) offers 
an alternative avenue to study the time evolution of the mass 
function. Figure [TT| shows the halo growth function for an 
8 Mpc box for three different starting redshifts, Zin = 50, 150, 
and 250 (these are the same simulations as in Fig. |6l). The re- 
sults are displayed at three redshifts, z = 20, 15, and 10 and 
for three mass bins, 10** - IO'^/i-'Mq, 10'^ - IO'^/j-'Mq, and 
lO^°-lO"/!-'M0. 

Assuming that the Warren fit scales at least approximately 
to high redshifts, the first halos in the lowest mass bin are pre- 
dicted to form at Zform ~ 25 (see Fig. |6j. We have found that 
if Zform is not sufficiently far removed from Zin, formation of 
the first halos is significantly delayed/suppressed. In turn, this 
leads to suppressions of the halo growth function and the mass 
function at high redshifts. As shown in FigurefTT] the suppres- 
sion can be quite severe at high redshifts: the simulation result 
at z = 20 from the late start at Zin = 50 is an order of magnitude 



lower than that from Zm = 250. At lower redshifts, the discrep- 
ancy decreases, and results from late-start simulations begin to 
catch up with the results from earlier starts. Coincidentally, the 
suppression due to the late start at Zin = 50 is rather close to the 
PS prediction which is very significantly below the Warren fit 
in the mass and redshift range of interest (see Fig.fTTTl. We take 
up this point further below. 




Redshift 



Fig . 11. — Halo growth function for an 8 h Mpc box stalled from three dif- 
ferent redshifts. The blue data points results from the z = 50 start, the turquoise 
data points from the z = 1 50 start, and orange from the z = 250 start, which is 
the redshift satisfying our starting ciiteria. The two fits shown are the War- 
ren fit (solid line) and the PS fit (dashed line). Three different mass bins are 
shown. It is interesting to note that the late start seems to follow the PS fit at 
high redshift. 



5.6. Finite-Volume Corrections 

The finite size of simulation boxes can compromise results 
for the mass function in multiple ways. It is important to keep 
in mind that finite-volume boxes cannot be run to lower than 
some redshift, zgnai, the stopping point being determined by 
when nonlinear scales approach close enough to the box size. 
Approaching too near this point delays the ride-up of nonlinear 
power towards the low-A; end, with a possible suppression of the 
mass function. 

As a consequence of this delay, the evolution (incorrectly) 
appears more linear at large scales than it actually should, as 
compared to the P{k) obtained in a much bigger box. There- 
fore, verifying linear evolution of the lowest A;-mode is by itself 
not sufficient to establish that the box volume chosen was suf- 
ficiently large. For all of our overlapping-volume simulations 
we have checked that the power spectra were consistent across 
boxes up to the lowest redshift from which results have been 
reported (Table 1 lists the stopping redshifts). 

Aside from testing for numerical convergence, it is important 
to show that finite-volume effects are also under control, es- 
pecially any suppression of the mass function with decreasing 
box size (due to lack of large-scale power on scales greater than 
the box size). Several heuristic analyses of this effect have ap- 
peared in the literature. Rather than rely solely on the unknown 
accuracy of these results, however, here we also numerically 
investigate possible systematic differences in the mass function 
with box size. 

Over the redshifts and mass ranges probed in each of our sim- 
ulation boxes, we find no direct evidence for an error caused by 
finite volume (at more than the ~ 20% level), as already em- 
phasized previously in iHeitmann et al. 1 (l2006al) . (Overlapping 
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Fig. 10. — Mass function at five different redshifts (z = 0, 5, 10, 15, and 20; top to bottom) compared to different fitting formulae. Note that the mass ranges are 
different at different redshifts. The simulation results have been corrected for FOF bias following Warren but not for finite-volume effects (for these, see Fig. 12). 
The bottom panel shows the ratio with respect to the Warren fit. Our simulations agree with the Warren fit at the 10% level for redshifts smaller than 10, although 
there is a systematic offset of 5% at z = 0, where our numerical results are higher than the fit. At higher redshifts, the agreement is still very good (at the 20% level) 
and becomes very close once finite- volume corrections are applied (Fig. ll2t . PS is a bad fit at all redshifts, and especially at high redshifts, where the difference 
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b ox-size results over diffe rent mass ranges are shown in Fig. 3 
of lHeitmann et al.]|2006ai ) Figure [TO] shows the corresponding 
results in the present work. This is not to say that there are 
no finite-volume effects (the very high-mass tail in a given box 
must be biased low simply from sampling considerations) but 
that their relative amplitude is small. Below we discuss how to 
correct the mass function for finite box size. 

5.6.1. Volume Corrections from Universality 

Let us first assume that mass function universality holds 
strictly, in other words, that for any initial condition the num- 
ber of halos can be described by a certain scaled mass func- 
tion (eq. HI) in which (j{M) is the variance of the top-hat- 
smoothed linear density field. In the case of infinite simula- 
tion volume, (j{M) is determined by equation (O, and the mass 
function F(M) of equation (|2|i is 



F(M): 



dn 



dlogM 



Ph t/lno- ' 

— /(cr) 

M-^ VlogM 



(46) 



In an ensemble of finite-volume boxes, however, one necessar- 
ily measures a different quantity: 



F'{M') 



dn' Ph I dXvLo' ' 



_ = ^/(a')l 

V Ml-' ' , 



(47) 



t/logM' M'' dlogM' 

Here a'{M') is determined by the (discrete) power spectrum of 
the simulation ensemble, although if universality holds as as- 
sumed, / in equations ( |46] l and WT\ is the same function. 

Since we are, in general, interested in the mass function 
which corresponds to an infinite volume, we can then correct 
the data obtained from our simulations as follows: for each box 
size we can define a function M'{M) such that 

(t(M) = (t'(M'(M)). (48) 

Using equations ( |46] ) - ( |48] |. we determine F(M) as 



F(M) = F'(M')—^. 

dM 



(49) 



Thus, the corrected number of halos in each bin is calculated as 



dn = dn' — 
M 



(50) 



The universality must eventually break down for sufficiently 
small boxes or high accuracy because the nonlinear coupling of 
modes is more complicated than that described by the smoothed 
variance. This violation can be partly corrected for by modify- 
ing the functional form of a'(M'). Therefore, we also explore 
other choices of a'(M') which may better represent the mass 
function in the box. To address this question we provide a short 
summary of the Press-Schechter approach. 

5.6.2. Motivation from Isotropic Collapse 

We first consider the idealized case of a random isotropic 
perturbation of pressureless matter and assume that the primor- 
dial overdensity at the center of this perturbation has a Gaus- 
sian probability distribution. The probability of local matter 
collapse at the center is then fully determined by the local 
variance of the primordial overdensity a^. Consequently, for 
the isotropic case the contribution of Fourier modes of vari- 
ous scales to the collapse probabihty is fully quantified by their 
contribution to a^. 

To see this, consider the evolution of matter density pioc at 
the center of the spherically symmetric density perturbation. 



For transparency of argument, let us focus on the evolution dur- 
ing the matter-dominated era; it is straightforward to general- 
ize the argument to include a dark energy component /9de(z), 
homogeneous on the length scales of interest, by a substitution 
Pioc Pm.ioc+/Ode in cquations ( fSTT i and ( |53] ). By Birkhoff 's law, 
the evolution of pioc and the central Hubble flow H\oc = 5^ ' ^loc 
are governed by the closed set of the Friedmann and conserva- 
tion equations, 

i ^ 



Fl]^^ — 



(51) 



a\o 



Po 

Ploc 



1/3 



^loc 

da 



loc 



dt 



(52) 



where k is a constant determined by the initial conditions, po is 
arbitrary (e.g., po = pz,|_o), and f is the proper time. 

The degree of nonlinear collapse at the center can be quanti- 
fied by a dimensionless parameter 



1- 



(53) 



First consider early times, when the evolution is linear, and let 
P\oc = Pbi^ Then for the growing perturbation modes dur- 
ing matter domination H\oc = H{1 - 6/3). Given these initial 
conditions, which set the initial pioc and the constant k in equa- 
tion ( ISTl i. the subsequent evolutions of pioc, i/ioc, and therefore 
q are determined unambiguously. 

During the linear evolution in the matter era q = 56/3 is small 
and grows proportionally to the cosmological scale factor a. For 
positive overdensity, nonlinear collapse begins when q becomes 
of order unity, reaching its maximal value q = I when H\oc = 0, 
and decreasing rapidly afterwards. (We can observe the latter 
by rewriting eq. 1531 as 

3k 



1/3 



(54) 



having applied eqs. lISTI and ll52l .) Nonlinear collapse of matter 
at the center of the considered region can be said to occur either 
when ^ ^ or when q reaches a critical "virialization" value qc. 

Now it is easy to argue that in the isotropic case the Press- 
Schechter approach gives the true probability of the collapse, 
P{q > qc,z), for a redshift z- Indeed, the evolution of q is set 
deterministically by the primordial density perturbation at the 
center; for adiabatic initial conditions specifically, it is set by 
the curvature perturbation ( at the center. Since higher values 
of ( lead to earlier collapse. 



P{q>qc,z)-- 



:P(C>C(z))=^erfc 



y/2a. 



(55) 



where the last equality uses the explicit form of P(C/) as a Gaus- 
sian distribution with a variance a^. 

If the considered isotropic distribution is confined by a (spher- 
ical) boundary and a at the center is reduced by removal of 
large-scale power, then equation ( |55] ) should accurately de- 
scribe the corresponding change of the collapse probability. In 
numerical simulations, due to the imposition of periodic bound- 
ary conditions, there is no power on scales larger than the box 
size. In this case the variance a should be specified by the ana- 
logue of equation (O with the integral replaced by a sum over 
discrete modes. 

For the mass function (eq. 1461 ), a constant reduction of the 
variance a^{M) due to the removal of large-scale power leads 
to a suppression of the mass function at the high-mass end and. 
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Fig. 12. — Mass function data corre cted fo r finite box volume by the ex- 
tended Press-Schecliter prescription of j|5.6.3l (squares). We show the res ults 
as a ratio with respect to the Warren fit and follow the conventions of Fig. 1101 
We also display the volume-uncon'ected data (crosses). Note that the volume- 
corrected data join smoothly across the box-size boundaries. This box correc- 
tion brings the results very close to universal behavior at high redshifts (see 

FigQD. 



counterintuitively, a boost at the low-mass end. The latter is 
easily understood as follows: The (T-dependent terms of equa- 
tion ( |46] |, 

^^"^^1^ = ^1^' ^^^^ 

give the fraction of the total matter density that belongs to the 
halos of mass M. When the variance is decreased by the box 
boundaries, this fraction is boosted at low masses due to a shift 
of halo formation to an earlier stage, where a larger fraction of 
matter is bound into low-mass objects. 



5.6.3. Numerical Results and Comparisons 

Following the abov e intuition, we emp loy the extended Press- 
Schechter formalism (iBond et al. Il99lh to correct for the miss- 
ing fluctuation variance on box scales. This formalism, while 
clearly inadequate at vari ous levels in describing halo formation 
in realistic simulations dBond et al. I Il99ll; iKatz et al. 1 119931: 
IWhite ..1996.) . has nevertheless been very successful as a cen- 
tral engine in describing the s tatistics of co smological struc- 
ture formation. As shown by iMo & White I (il996) using A^- 
body simulations, the biasing of halos in a spherical region 
with respect to the average mass overdensity in that region is 
very well described by th e extended Press-Schechter approach. 
iBarkana & Loeb I (|2004|) discussed the suppression of the halo 
mass function in terms of this bias, and suggested a prescription 
for adjusting large-volume mass function fits such as Warren or 
ST to small boxes. Here we do not follow this path but directly 
work with the numerical data by correcting the number of halos 
in each bin as in equation (ISOl l. 

In the extended Press-Schechter scenario of halo formation, 
a' on the right-hand side of equation WT\ would be approx- 
imately connected with a via a'^ = - dBond et al. I 

where o\f^^^ is the variance of fluctuations in spheres 
that contain the simulation volume. Since extended Press- 
Schechter theory is derived for spherical regions, while our sim- 
ulation boxes are cubes, we define /?(box) as the radius of a 
sphere enclosing the same volume as in the simulations. 

The action of this correction is shown in Figure [121 Finite- 
volume corrections are subdominant to statistical error at z = 
and 5. At higher redshifts, the corrections produce results 
that are consistent across box sizes, i.e., that have no systematic 
shape changes or "jumps" across box boundaries. Moreover, 
the action of the corrections is to bring the simulation results 
closer to a universal behavior. We discuss this aspect further 
below. 

For completeness, we mention two other approaches aimed 
at b ox-ad justing the m ass function. The first ( Yoshida et al. I 
l2003cl : rBagla & PrasadT2006h simply replaces the original mass 
variance (eq. (O) with 



1 



dHz) 



k^P(k)W^(k,M)dk 



(57) 



2ir/L 



the lower cut-off arising from imposing periodic boundary con- 
ditions (L is the box-size). (For enhanced fidelity with simula- 
tions, the integral in eq. ISTl goes to a sum over the simulation 
box modes.) This approach basically assumes that cr defined 
via an infrared cutoff is the appropriate replacement for the 
infinite-volume mass variance. Figure [13] shows the effect of 
this suggested correction: At z = and 5 it is not noticeable, but 
at higher redshifts the correction is significant relative to the 
accuracy with which the binned mass function is determined. 
Furthermore, it exhibits systematic shape changes and offsets 
across boxes, in contrast to the results shown in Figure [12] For 
example, at z = 10 the corrected data at the crossover point be- 
tween the 4 and 8/!"'Mpc boxes (~ IO'^H'^Mq) have an offset 
of 5%. We conclude that this approach is disfavored by our 
simulation results. 

An alternative strategy is to estimate the mass variance from 
each realization of P(k) in the individual s imulation boxes an d 
to treat every box individually, as done in lReed et al. I ( l2007h . 
This has in fact two purposes: to compensate for the reaUzation- 
to-realization variation in density fluctuations (which could be a 
problem for small boxes) and also to compensate for an overall 
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Fig. 13.— 

Mass function corrected for a finite box using the assumption of 
strict universality, as described in ^5.6.1| (squares). Again, we 
show uncorrected data as well (crosses), and follow the con- 
ventions of Fig. [To] This correction produces a clear systematic 
shift in the results across box boundaries. 



suppression in the mass function as discussed above. The dis- 
advantage is that each of many realizations now has a different 
cr(M) for a given value of M. 

5.7. Mass Function Universality 

Finally, we investigate the universality of the mass func- 
tion found by Jenkins. Approximate universality is expected 
from the analytic arguments of PS and the extended, excursion- 
set formulation of Bond et al. (1991). The universal be- 
havior of halo formation persists even in the model of ellip- 
soidal collapse of ST, in which the predicted mass function 



is no longer of the PS form. On the other hand, the univer- 
sality cannot be exact if the nonlinear interactions of differ- 
ent scales are fully accounted for: The nonlinear evolution 
that leads to the formation of halos of mass M must involve 
multiple degrees of freedom that are described by more pa- 
rameters than the overall variance of the primordial overden- 
sity smoothed by a top-hat filter W(r,M). The universality is 
expected to be violated at sufficiently high resolution of the 
mass function even in the PS-type spherical collapse model: 
It is more reasonable to represent the probability of the col- 
lapse not by a fraction of particles at the center of spheres en- 
closing a mas s M but by any fraction of particles belong ing to 
such spheres dBetancort-Riio & Montero-Dorta I l2006al) . The 
improved mass-function derive d from this argument deviates 
somew hat from a universal form dBetancort-Riio & Montero-Dorta I 

HoolB). 

To investigate the extent our numerical simulations are con- 
sistent with universality, we combine our results for /((T,z) as 
a function of the variance cr~' from the entire simulation set in 
one single curve at various redshifts. This curve is expected 
to be independent of redshift if universality holds. We display 
the results in Figure [14] for the raw data and in Figure [15] for 
the same data after applying the volume corrections discussed 
earlier 

In the raw data of Figure [14] the agreement with the various 
fits is quite tight (except for PS) until In cr~' > 0.3. Beyond this 
point, the multiple-redshift simulation results do not lie on top 
of each other; in the absence of any possible systematic devia- 
tion, this would denote a failure of the universality of the FOF, 
h = 0.2 mass function at small cr. Note also that beyond this 
point the ST and Jenkins fits have a steeply rising asymptotic 
behavior (relative to the Warren fit). The Reed et al. (2003) fit, 
meant to be valid over the range -1 .7 < In cr"' < 0.9, is in better 
agreement with our results, to the extent that a single fit can be 
overlaid on the data. 

The ostensible violation of universality seen above is small, 
however, and subject to a systematic correction due to the finite 
simulation volume(s). On applying the correction discussed in 
Section 5.6.3, we obtain the results shown in Figure[T5] the key 
difference being that beyond Incr"' > 0.3 the multiple-redshift 
simulation results now lie on top of each other and, within the 
statistical resolution of our simulations, are consistent with uni- 
versal behavior Spe cifically, we do no t observe the sort of vi- 
olation reported by Ree d et al~l (|2007|) at high redshifts. This 
could be due to several factors. The finite-sampling FOF mass 
correction and the fi nite-volume correcti ons we employ are dif- 
ferent from those of iReed et al. 1 (I2007h and the boxes we use 
at high redshifts are significantly larger. We note also that the 
difference between the Warren fit and the z-dependent fit of 
iReed et al. I (|2007|) does not appear to be statistically very sig- 
nificant given their data. 

6. CONCLUSIONS AND DISCUSSION 

We have investigated the halo mass function from A^-body 
simulations over a large mass and redshift range. A suite of 60 
overlapping-volume simulations with box sizes ranging from 4 
to 256/!~'Mpc allowed us to cover the halo mass range from 
10^ to 1O'^^/!~'M0 and an effective redshift range from z = Q 
to 20. 

In order to reconcile conflicting results for the mass func- 
tion at high redshifts, as well as to investigate the reahty of the 
breakdown of the universality of the mass function, we have 
studied various sources of error in A^-body computations of the 




Fig. 14. — Scaled differential mass function from all simulations, prior to ap- 
plying finite-volume corrections. Fits shown are Wan'en (red), PS (dark blue), 
ST (black), Jenkins (light blue), and Reed et al. (2003) (yellow). Dashed lines 
denote an extrapolation beyond the original fitting range. The bottom panel 
shows the ratio relative to the Warren fit. The failure of the different redshift 
results to lie on top of each other at small values of a indicate a possible vio- 
lation of universality. 



mass function. A set of error control criteria need to be satis- 
fied in order to obtain accurate mass functions. These simple 
criteria include an estimate for the necessary starting redshift, 
for the required mass and force resolution to resolve the halos 
of interest at a certain mass and redshift, and for the number of 
time steps. 

The criteria for the initial redshift appear to be particularly 
restrictive. For small boxes, commonly used in the study of 
the formation of the first objects in the Universe, significantly 
higher initial redshifts are required than is the normal practice. 
A violation of this criterion leads to a strong suppression of the 
mass function, most severe at high redshifts. Recent results by 
other groups may be contaminated due to a violation of this 
requirement; a careful re-analysis of small-box simulations is 
apparently indicated. 

The force resolution criterion is especially useful for grid 
codes, PM as well as adaptive mesh. The mass function can be 
obtained reliably from PM codes down to small-mass and up 
to high-mass halos provided the halos are adequately resolved. 
The resolution criterion is also very useful in setting refinement 
levels for adaptive mesh refinement (AMR) codes. As shown 
recently (O'Shea et al. 2005; Heitmann et al. 2005, 2007) the 
mass function from AMR codes is suppressed at the low-mass 
end if the base refinement level is too coarse. A more detailed 
analysis on how to incorporate our results to improve the effi- 
ciency of AMR codes is underway. 

The results for the required number of time steps to resolve 
the mass functions is somewhat surprising. The halo mass func- 
tion appears to be very robust with respect to the number of 
time steps chosen to follow the evolution, even though the in- 
ner structure of the halos will certainly not be correct. Even 
a small number of time steps is sufficient to obtain a close-to- 
correct mass function at z = 0. This considerably simplifies the 
study of the mass function and its evolution. 



F ig. 15. — Volume-corrected scaled differential mass function following 
Fig[T4] Note the significantly improved agreement with universal behavior 
(overlapping results beyond Incr"' ~ 0.3). 



Since finite-volume effects can also lead to a suppression of 
the mass function, we have tried to minimize the importance of 
these effects by avoiding too-small box sizes, by using overlap- 
ping boxes, and by restricting the mass range investigated in a 
given box size. In addition, we have found that a box-size cor- 
rection motivated by the extended Press-Schechter formalism 
for the mass variance appears to give consistent results when 
applied to our multiple-box simulation ensembles. 

We now briefly comment on results found previously by 
other groups. Jang-Condell & Hernquist (2001) find good 
agreement with the PS fit at z = 10 for a mass range 4 x 10^ - 
4 X 10**/i~'Mq. The crossover of PS with the more accurate fits 
at z = 10 takes place in exactly this region (see Figs.fTlandfTOli. 
Therefore, all fits are very close, and the mass function from a 
single 1 Mpc box at a single redshift as shown in Jang-Condell 
& Hernquist (2001) cannot distinguish between them. 

As mentioned earlier in ^2.51 good agreement with the PS 
result has been reported at high redshifts (some results being 
even lower than PS) by seve ral other groups (Yoshida et al. 
2003a. 2003b. 2003c;ICen et al. 2004; Trac &_Cer i 2006). The 
simulations of ICen et al. I ( |2004|) and ffrac & Cen I ilOOH) were 
started at Zin ~ 50, substantially below the starting redshift that 
would be sugg ested by our work. The very large number of 
particles in the Trac & Cen I ([2006) simulation requires a high 
starting redshift (Fig.|5]l. Therefore, the depressed mass func- 
tion results of these s imulations are very c onsistent with a too- 
low initial redshift. dTrac & Cen I (1200 6') have recently rerun 
their simulations with a much higher initial redshift [z = 300], 
and now find results con sistent with ours.) The initial particle 
density of the llliev et al. 1(1200 6) simulations is very close to that 
of our 16/z~'Mpc box, in which case also a high redshift start is 
indicated (we used Zm = 200). Finally, the initial redshift of the 
Yoshida et al. papers, Zin = 100, for boxes of size ^ 1 /i"'Mpc, 
also appears to be significantly on the low side. 

We have compared our simulation results for the mass func- 
tion with various fitting functions commonly used in the liter- 
ature. The recently introduced (z = 0) fit of Warren leads to 
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good agreement (at the 20% level with no volume correction, 
and at the 5% level with volume correction) at all masses and 
all redshifts we considered. Other modern fits, such as Reed et 
al. (2003, 2007), also lie within this range. These fits do not 
suffer from the overprediction of large halos at high redshifts 
observed for the ST fit. The PS fit performs poorly over almost 
all the considered mass and redshift ranges, at certain points 
falling below the simulations by as much as an order of magni- 
tude. 

The evolution of the mass function can be used to test the 
(approximate) universality of the FOF, b = 0.2 mass function. 
At low redshifts o ur data are in good agreement with those of 
iReed et al. I (12007 ) (at z = 5), finding a (possible) mild redshift 
dependence (at the 10% level). At higher redshifts, however, 
we find that volume corrections are important to the extent that 
little statistically significant evidence for breakdown of univer- 
sality remains in our mass function data. A full theoretical un- 
derstanding of this very interesting result remains to be eluci- 
dated. 

We have made no attempt to provide a fitting function for 
our data due to several reasons. First, the current simulation 
state of the art has not reached the point that one can be confi- 
dent of percent-level agreement between results from different 
simulations even in re gimes that are not statistics-dominated 
dHeitmann et al. I l2007h . Second, simulations have not suffi- 
ciently explored the extent to which universal forms for the 
mass function are indeed applicable as cosmological parame- 
ters are systematically varied. Third, absent even a compelling 
phenomenological motivation for the choice of fitting func- 
tions, there is an inherent arbitrariness in the entire procedure. 
Finally, it is not clear how to connect the FOF mass function 
to observations. In general, tying together mass-observable 
relations requires close coupling of simulations and observa- 
tional strategies. In studies of cosmological parameter estima- 
tion, we support working directly with simulations rather than 
with derived quantities, which would add another layer of pos- 
sible systematic error. Because observations already signifi- 
cantly constrain the parametric range, and are a smooth func- 
tion of th e parameters, this approach is quite viable in practice 
dHeitmann et a l. 2006b: Habibeta l. 200% . 
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APPENDIX 

In this Appendix we discuss in detail previous results on the 
mass function at high redshift. As explained in the main paper, 
these results are often contradictory. We structure our discus- 
sion with respect to the physical volume simulated. 

Small-Volume Simulations 

Small-box simulations of side ^ l/!"'Mpc have been per- 
formed by several groups. Using a treecode with softening 
length 0.4/i~'kpc, and a l /i~'Mpc box with 128^ particles, 
iJangj-Condell & Hernquist i 1200 1 ) evolved their simulation from 
Zin = 100 to z = 10. With a halo finder that combined overdensity 
criteria with an FOF algorithm, the mass function was deter- 
mined over the range 10^^- 10^ ' H'^Mq, keeping halos with as 
few as eight particles. At z = 10 they found "remarkably close 
agreement" with the PS fit but did not quantify the agreement 
explicitly. 

In a series of papers, Yoshida et al. ran simulations with sim- 
ilar box sizes as above, most including the effects of gas dynam- 
ics. The simulations were performed with the TreePM/smoothed 
particle hydrodynamics code GADGET-II ( Springel 2005) and 
followed the evolution of 2 x 324^ particles (324^ in the case of 
dark matter only), covering a halo mass range of IO^-IO'^Mq. 
All simulations were started at Zm = 100 from "glass" ini- 
tial conditions (Baugh et al. 1995; White 1996), in contrast 
to the grid-based init ial conditions used here. The focus of 
lYoshida et al. I (l2003al) was the origin of primordial star-forming 
clouds. As part of that investigation, a dark-matter-only simu- 
lation in a 1.6/;~'Mpc box was carried out. The halo density 
results for z = 20 to 32 lay systematically below the PS predic- 
tion, with the discrepancy being worse at high redshifts. The 
authors argued that this low abund ance of halos was (poss ibly) 
due to finite-box-size effects. In Yoshida et al. I d200 3b'). the 
mass function at z = 20 for a warm dark matter model was com- 
pared with C DM, with the s i mulatio n set up being very sim- 
ilar to that of lYoshida etal.l d2003al) . a 1 Mpc box started at 
z = 100. The results obtained were also similar; at z = 20 the 
CDM mass function was in g ood agreement with the PS fit. In 
a third paper, Yoshida et al. I a2003 c), a running spectral index 
was considered. Here results for a standard CDM mass function 
for a 1 Mpc box were given, this time at z = 17 and 22. Con- 
sistent with their previous results, they found good agreement 
with PS at these redshifts. (The FOF linking length used in the 
last paper was b = 0.2, while in the first two papers b = 0.164 
was chosen. This did not appear to make much of a difference, 
however.) These papers do not quantitatively compare the nu- 
merical mass function to the PS fit. (In contrast to these find- 
ings, a recent 1 Mpc b ox GADGET-II sim ulation with Zin ^120 
has been performed bv lMaio et al. 1 (12006) who find good agree- 
ment with the Warren fit as extrapolated by linear theory - in 
clear disagreement with PS.) 

A similar strategy was followed in ICen et alT 1 (120041) who 
investigated dark matter halos in a mass range of 10^'^ to 
lO'^r'Mo, using a TreePM code (Xu 1995; Bode et al. 2000). 
The box size was taken to be 4/;"' Mpc, the softening length was 
set at 0.14/!"'kpc, 512-' particles were used, and the simulations 
had a starting redshift of Zin = 53. Hal os were identified us- 
ing the overdensity scheme DENMAX dBertschinger & Gelb I 
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Il99lh . Among other quantities, they studied the mass function 
between z = 1 1 and 6 and found that the PS function "provides 
a good fit" but without expHcit quantification. 

Overall, these small-box simulations, run with different codes 
and different halo finders, all found a "depressed" mass func- 
tion (see Fig.[T]), consistent with PS and deviating very signif- 
icantly from the predictions of the more modern fitting forms. 
In contrast, other simulations also using small boxes hav e come 
to qui te different conclusions. For example, in Reed et al. I 
(l2007h . a large suite of different box sizes and simulations was 
used to cover the mass range between 10^ and 1O"'^/!~^M0 at 
high redshift. The smallest boxes considered in this study were 
1 /z"'Mpc on a side. The authors studied the halo mass function 
at redshifts out to z = 30, implementing a correction scheme to 
account for finite-box effects, a s discussed in more detail be- 
low. Overall, Reed et al. (2007) confirmed p revious resu lts as 
found by Reed et al. (2003) and Heitmann etd"] (|2006a|) : PS 
underestimates the mass function considerably (by at least a 
factor of 5 at high redshift and high masses), and ST overpre- 
dicts the halo abundance at high redshift. 

Large-Volume Simulations 

The large-box strategy is exemplified by a recent dark matte r 
simulation with the GADGET-II code (" Springel etani2005h . 
The evolution of 2160^ particles in a 500/j-^Mpc box was fol- 
lowed from Zin =127 until z = 0. The softening length was 
5/i~'kpc. The high mass and force resolution was sufficient 
to study the mass function reliably down to a redshift of z = 10, 
covering a mass range of 10"^ to 10"'/!~'Mq, with halos being 
identified by a standard FOP algorithm with b = 0.2. The results 
are consistent with the Jenkins fit, even though the mass func- 
tion points at redshifts z = 1 .5, 3.06, and 5.72 are slightly higher 
than the Jenkins fit and slightly lower for z = 10. No residuals 
were shown nor quanti tative statements m ade. 

In two recent papers. llliev et al. I (l2006h and lZahn et"ain ( l2007h 
investigated cosmic reionization, providing mass funct i on re- 
sults at high redshift as part of this work. Illiev et aL 1 (120061) 
ran a PM simulation with PMFAST (Merz et al. 2005) in a 
100/z"'Mpc box with 1624^ particles on a 3248^ mesh. They 
present results for the mass function at redshifts between z = 6 
and 18.5, using a spherical overdensity halo finder. At lower 
redshifts they find good agreement with ST, and at high redshift 
(z > 10) the results are closer to PS (because of their limited 
mass range, a mor e quantitative statement is difficult to make). 
IZahn et al. I (l2007 t) ran a 1024"' particle simulation (dark matter 
only) in a 65.6/z~'Mpc box with GADGET-II and analyzed the 
FOF, b = 0.2 mass function out to z = 20. Between z = 6 and 
14 they found good agreement with ST in the mass range of 
10^ to IO^^Mq. At z = 20 they found that the simulation results 
were below ST but above PS, in relatively good agreement with 
the rec ent findings of .Heitmann et al. i (i2006a) and iReed et al. i 
(l2007h . 

Medium Volume Simulations 

iReed et aTl (l2003h chose a compromise between the large- 
and small-box strategies by picking a 50/i~'Mpc box sampled 
with 432^ particles. The tree code PKDGRAV was used to 
evolve the simulation from different starting redshifts between 
Zin = 139 and 69 until z = 0. The smallest halo contained 75 par- 
ticles, leading to a mass range of roughly 10'" to 10'"*^/!~'Mq. 
Good agreement (better than 10%) was found with the ST fit 
up to z — 10. For higher redshifts, the ST fit overpredicted the 



number of halos, up to 50% at z = 15. At this high redshift, 
statistics were lacking, and the resolution was not sufficient to 
resolve very small halos. A more recent 50/i~'Mpc simulation 
with PMFAST with Zin = 60 has been carried out by Trac & CenI 
(2006) using a spherical overdensity definition of halo mass. In 
this work, the mass function, in the redshift range 6 < z < 15, is 
found to be in very good agreement with PS, in gross contradic- 
tion with the results of most of the other simulations mentioned 
above. (This contradiction has recently been resolved by rerun- 
ning their simulation with z,„ = 300 and identifying halos with 
a /7 = 0.2 FOF finder) 
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